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Abstract 

This  paper  reports  on  the  systematic  experimental  validation  of  a  comprehensive  3D  CFD-based  computational  model  presented  and  documented 
in  Part  1.  Simulations  for  unit  cells  with  straight  channels,  similar  to  the  Ballard  Mk902  hardware,  are  performed  and  analyzed  in  conjunction 
with  detailed  current  mapping  measurements  and  water  mass  distributions  in  the  membrane-electrode  assembly.  The  experiments  were  designed  to 
display  sensitivity  of  the  cell  over  a  range  of  operating  parameters  including  current  density,  humidification,  and  coolant  temperature,  making  the 
data  particularly  well  suited  for  systematic  validation.  Based  on  the  validation  and  analysis  of  the  predictions,  values  of  model  parameters,  including 
the  electro-osmotic  drag  coefficient,  capillary  diffusion  coefficient,  and  catalyst  specific  surface  area  are  determined  adjusted  to  fit  experimental 
data  of  current  density  and  ME  A  water  content.  The  predicted  net  water  flux  out  of  the  anode  (normalized  by  the  total  water  generated)  increases 
as  anode  humidification  water  flow  rate  is  increased,  in  agreement  with  experimental  results.  A  modification  of  the  constitutive  equation  for  the 
capillary  diffusivity  of  water  in  the  porous  electrodes  that  attempts  to  incorporate  the  experimentally  observed  immobile  (or  irreducible)  saturation 
yields  a  better  fit  of  the  predicted  MEA  water  mass  with  experimental  data.  The  specific  surface  area  parameter  used  in  the  catalyst  layer  model  is 
found  to  be  effective  in  tuning  the  simulations  to  predict  the  correct  cell  voltage  over  a  range  of  stoichiometries. 

©  2008  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

A  number  of  reviews  have  recently  provided  a  status  of  com¬ 
putational  fluid  dynamics  (CFD)  models  for  PEMFCs  [1-4]. 
One  of  the  main  focuses  in  this  paper  is  the  issue  of  vali¬ 
dation.  Table  1  summarizes  recent  modeling  and  simulation 
studies  on  PEMFC  and  the  methods  used  in  validating  them. 
Until  recently,  most  numerical  simulations  relied  on  validating 
computational  results  by  so-called  zero-dimensional  data  such 
as  the  V-I  curve  and  overall  water  balance,  which  are  global, 
spatially  averaged  values.  This  approach  provides  limited  con¬ 
fidence,  even  when  the  number  of  data  used  for  validation  is 
abundant;  for  instance  entirely  different  electric  field  distribu¬ 
tions  obtained  from  simulations  in  which  the  electrochemical 
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asymmetry  factor  is  varied  can  result  in  essentially  identical 
global  polarization  curves  [5] .  Studies  using  reduced  dimensions 
and  simplified  models  and  listed  in  Table  1,  focus  on  analysis 
(a)  across  the  MEA  [6],  (b)  along-the-channel  [7-13],  and  (c) 
under-the-rib  [14].  Each  approach  has  its  merits  in  providing 
physical  insights  into  the  transport  phenomena  in  the  geometric 
space  considered.  The  validation  with  global  data  such  as  polar¬ 
ization  curves  or  water  balance,  however,  should  be  interpreted 
with  caution  as  there  are  uncertainties  associated  with  some  of 
the  underlying  assumptions  in  the  simplified  models  as  well  the 
multi-dimensional  effects  not  considered  in  the  models.  When 
locally  resolved  data  are  used  for  validation  [15,16],  the  same 
concern  remains.  The  uncertainties  due  to  model  simplification 
and  reduced  dimensions  are  obviously  relaxed  when  compre¬ 
hensive,  3D  simulations  are  employed  [17-23],  but  uncertainties 
remain  due  to  fitting  of  global  data.  More  recently,  some  numer¬ 
ical  simulations  have  been  validated  against  local  experimental 
data  [24-26]. 


424 


PC.  Sui  et  al.  /  Journal  of  Power  Sources  180  (2008)  423^-32 


Nomenclature 

Cp  specific  heat  ( J  mol- 1  K- 1 ) 

Di  water  diffusion  coefficient  (mol  m- 1  s  - 1 ) 

I  current  density  (A  m-2) 

m  mass  flow  rate  of  fluid  (kg  s-1) 

rid  electro-osmotic  drag  coefficient 

P  pressure  (Pa) 

RH  relative  humidity 

.v  saturation 

S/V  specific  surface  area  (m2  m-3) 

T  temperature  (K) 

Greek  letter 
X  water  content 

Subscripts 
A  anode  side 

C  cathode  side 


In  the  present  study,  the  comprehensive  3D  CFD  code 
described  in  Part  1  [27]  is  validated  against  experimental  data 
obtained  for  a  unit  cell.  The  data  include  local  current  density 
distribution  and  mass  of  water  along  the  channel  measured  over  a 
wide  range  of  operating  conditions  that  results  in  significant  dif¬ 
ferences  in  the  dominant  processes  and  couplings,  making  this 
data  set  particularly  challenging  for  3D  simulations.  The  objec¬ 
tives  of  this  paper  are  two-fold:  first,  validation  of  the  physical 
models  and  some  determination  of  model  parameters  from  the 
experimental  data  set,  and,  second,  parametric  analysis  to  gain 
insight  into  the  extent  and  effects  of  salient  transport  processes. 

2.  Methodology 

2.7.  Experimental  data 

The  experimental  data  used  for  validation  were  collected 
using  the  MRED  (MEA  Resistance  and  Electrode  Diffusion) 
method  developed  by  Stumper  et  al.  [28].  This  method  allows  in 
situ  determination  of  MEA  resistance  and  electrode  diffusivity 
of  the  cell.  The  anode  and  cathode  surfaces  of  the  unit  cell  are 
attached  to  16  current  collector  pucks,  which  connect  to  the  elec¬ 
trical  load.  The  current  through  each  puck  is  determined  by  the 
voltage  drop  across  a  shunt  resistor  on  each  puck.  The  test  cell 
is  operated  on  a  custom-designed  test  station  allowing  accurate 


Fig.  1.  Comparison  of  water  balance  data  and  numerical  predictions. 


control  and  monitoring  of  all  operating  parameters.  Steady- state 
polarization  curves  are  obtained  under  constant  fuel  and  oxidant 
stoichiometry  with  respect  to  the  total  current. 

The  mass  of  water  in  the  MEA  is  determined  using  a  non¬ 
destructive  method  that  measures  the  weight  of  water  ex  situ. 
The  mass  of  the  water  is  measured  individually  for  each  of 
the  16  pucks  to  provide  the  water  distribution  for  the  MEA. 
Experimental  data  are  collected  for  a  range  of  operating  con¬ 
ditions.  These  data  are  used  to  validate  the  model  predictions 
by  post-processing  the  3D  CFD  predictions  in  exactly  the  same 
manner  as  in  the  experiments.  Post-processing  of  the  3D  fields 
is  described  in  Part  1  [27] . 

3.  Results  and  discussion 

3.1.  Water  balance:  numerical  versus  experimental 

Water  balance,  defined  as  the  net  water  transfer  between  the 
anode  inlet  and  outlet  normalized  by  the  total  water  generated  by 
the  cell,  is  a  measurable  quantity  that  characterizes  water  man¬ 
agement  of  a  unit  cell.  A  negative  value  indicates  a  net  transfer 
from  anode  to  cathode.  As  far  as  validation  of  numerical  simula¬ 
tions  is  concerned,  the  water  balance  is  a  zero-dimension  integral 
quantity,  and  should  therefore  be  interpreted  with  caution  regard¬ 
ing  local  water  transfer  across  the  membrane  inside  the  unit 
cell,  in  particular  for  a  counter-flow  configuration  when  water 
may  diffuse  from  cathode  to  anode.  Fig.  1  shows  a  comparison 


Table  1 


Summary  of  literature  on  numerical  simulation  and  validation 

Model 

Data  used  for  validation 

Global  data:  V-I,  water  balance 

Local  data:  I,  T,  y;,  mw 

Reduced  dimension, 
simplified  model 
Comprehensive,  3D 

Bernardi  and  Verbrugge  [6],  Gurau  and  Liu  [7],  Yi  and  Nguyen  [8],  Rowe  and  Li  [9],  You  and  Liu  [10], 
Siegel  et  al.  [11]  and  Kulikovsky  [15] 

Um  and  Wang  [16],  Dutta  et  al.  [17]a,  Zhou  and  Liu  [18],  Berning  and  Djilali  [19]a,  Lee  et  al.  [20]a, 
Mazumder  and  Cole  [21],  Li  and  Becker  [23]  and  Sivertsen  and  Djilali  [5] 

Berg  et  al.  [13]  and 
Kulikovsky  [15] 

Ju  and  Wang  [24]  and 
Li  and  Becker  [23] 

The  membrane  was  not  spatially  resolved  in  the  computation. 
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Table  2 


Summary  of  experimental  conditions  for  CFD  input 


Case 

Nominal  specs 

I  (A  cm-2) 

Cathode 

Anode 

DP, ox  (C) 

Tox,in  (C) 

RH 

stoich 

DP,h2  (C) 

Th2,in  (C) 

RH 

stoich 

1 

Baseline 

1.0 

62.0 

65.9 

0.839 

1.66 

54.4 

74.2 

0.411 

1.50 

2 

RH50% 

1.0 

52.3 

65.1 

0.549 

1.71 

38.4 

73.9 

0.185 

1.54 

3 

RH25% 

1.0 

37.2 

65.3 

0.252 

1.76 

56.7 

68.3 

0.590 

1.52 

4 

RH0% 

1.0 

0.0 

65.4 

0.026 

1.78 

57.0 

62.3 

0.782 

2.24 

5 

1  =  0.45 

0.45 

62.3 

65.3 

0.874 

1.66 

54.7 

74.9 

0.405 

1.50 

6 

1  =  0.1 

0.1 

62.7 

64.5 

0.924 

2.82 

54.5 

74.5 

0.408 

1.50 

7 

Stoich  =1.6 

1.0 

61.9 

65.5 

0.849 

1.48 

54.1 

75.0 

0.390 

1.51 

8 

Stoich  =1.4 

1.0 

61.7 

64.7 

0.873 

1.29 

31.9 

75.4 

0.122 

1.57 

9 

Stoich  =1.2 

1.0 

61.9 

65.1 

0.865 

1.11 

0.0 

75.3 

0.017 

1.59 

Pressure  at  inlet:  Pa  =  3.2  bar,  Pc  =  3.0 bar. 


of  water  balance  versus  anode  humidification  water  flow  rate 
for  the  baseline  case  with  different  anode  inlet  dew  point  tem¬ 
peratures.  The  average  current  density  is  1  Am-2  for  all  cases 
shown.  Both  experiments  and  predictions  show  a  linear  varia¬ 
tion  of  water  balance  with  anode  humidification  water  flow  rate. 
Under  dryer  anode  conditions,  the  experimental  data  shows  neg¬ 
ative  water  balance,  suggesting  back  diffusion  from  the  cathode 
side. 

3.2.  Compilation  of  MEA  water  mass  data 

Experimental  data  obtained  by  the  MRED  method  for  a 
unit  cell  are  used  in  the  present  investigation.  The  test  con¬ 
ditions  are  summarized  in  Table  2.  The  MRED  data  basically 
consist  of  water  content  profiles  and  current  mapping  data  for 
unit  cell  operated  under  various  operating  conditions.  To  illus¬ 
trate  the  broad  range  of  operating  conditions  over  which  the 
data  was  taken,  Fig.  2  shows  all  MEA  water  mass  data  plotted 
versus  current  density  for  the  13  test  cases  listed  in  Table  2. 
Two  observations  from  this  figure  are:  (a)  water  content  of  the 
MEA  populates  the  range  of  5  zb  0.5  mg  cm-2 — this  is  consid¬ 
ered  as  a  “saturated  MEA”,  corresponding  to  a  fully  humidified 
membrane  and  a  partially  saturated  GDL  having  a  maximum 


saturation  of  ca.  0.2,  cf.  Appendix  A  and  (b)  low  current  den¬ 
sity  results  correlate  to  low  water  content,  suggesting  that  low 
current  density  conditions  might  be  associated  with  relatively 
higher  ohmic  resistance  in  the  MEA.  The  broad  range  of  oper¬ 
ating  conditions  documented  experimentally  make  this  data  set 
particularly  useful  in  assessing  the  performance  of  the  model 
over  a  large  operating  envelope  with  different  dominant  and/or 
limiting  transport  mechanisms. 

3.2.1.  Sensitivity  to  inlet  humidification 

Fig.  3  shows  a  comparison  of  experimental  data  versus  CFD 
baseline  results  for  four  humidification  cases,  i.e.  RH  for  84, 
55,  25  and  0.03%,  respectively.  One  can  see  that  the  numerical 
results  are  consistent  with  experimental  data  qualitatively  for 
different  humidification  but  there  are  quantitative  discrepancies 
between  them.  The  numerical  predictions  fail  to  show  the  sat¬ 
urated  MEA  region  in  the  middle  of  the  cell.  The  drop  off  of 
water  content  in  the  cathode  outlet  (high  puck  number)  appears 
to  be  more  pronounced  in  the  data  than  numerical  prediction. 

Several  attempts  of  adjusting  the  MEA  properties  are  taken 
to  fit  the  data:  (1)  capillary  diffusion  coefficient  used  for  porous 
media,  (2)  membrane  properties  including  the  sorption  isotherm, 
drag  coefficient  and  water  diffusivity,  and  (3)  specific  surface 
area  in  the  cathodic  catalyst  layer.  Comparison  and  analysis  of 
numerical  prediction  versus  experimental  data  are  given  in  the 
following  sections. 


Fig.  3.  MEA  water  content  profiles:  experimental  data  versus  CFD  baseline. 
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3.3.  Capillary  diffusion  coefficient 

The  rationale  in  changing  the  capillary  diffusion  coefficient 
is  that  a  number  of  porous  media  exhibit  an  immobile  or  crit¬ 
ical  saturation  value  (also  referred  to  as  irreducible  saturation) 
below  which  the  effective  capillary  diffusion  is  zero.  Very  recent 
pore  network  simulations  targeted  at  GDL  media  exhibit  criti¬ 
cal  saturation  values  in  the  range  of  0. 1-0.2  [29,30].  This  is 
consistent  with  the  inference  from  the  results  that  the  default 
capillary  diffusion  coefficient  appears  to  high,  thus  resulting  in 
lower  than  observed  water  retention  in  the  GDL.  A  numeri¬ 
cal  experiment  was  thus  performed  to  investigate  the  effects  of 
capillary  diffusion  coefficient  on  predicted  GDL  water  content. 

Fig.  4  shows  a  comparison  of  the  default  and  modified  cor¬ 
relations  for  the  capillary  diffusion  coefficient.  The  modified 
correlation  tested  here  has  a  trend  similar  to  the  default  and  is 
represented  by  two  segments  over  the  range  v  =  [0,1],  but  the 
overall  values  are  lower  than  in  default  correlation.  The  formula 
for  the  tested  case  takes  an  exponential  form: 

£>cap  =  A  •  [exp (ks)  -  1]  (1) 

For  s < 0.2,  A  =  10-5,  £  =  57.6;  for  s>0.2,  A=  1(T2,  £=5.77. 
The  low  saturation  region  in  the  tested  correlation  has  a  region 
of  very  low  value  to  mimic  the  immobile  regime  for  water. 
Fig.  5  shows  the  predicted  mass  of  water  using  the  hypothet¬ 
ical  correlation  of  capillary  diffusion  coefficient  for  the  same 
cases  presented  in  Fig.  3.  One  can  see  that  the  mass  of  water  in 
the  middle  portion  is  shifted  up  to  be  closer  to  the  experimental 
data,  while  on  both  ends  of  the  cell  the  discrepancies  remain. 
It  is  noted  that  the  boundary  dividing  the  single-phase  regime 
and  two-phase  regime  is  not  affected  by  the  correlation  used  for 
capillary  diffusion  coefficient. 

The  effects  of  capillary  diffusion  coefficient  can  be  seen  in 
Fig.  6,  which  shows  a  more  diffuse  distribution  of  saturation  in 
the  GDL  using  the  default  correlation.  With  reduced  value  of 
capillary  diffusion  coefficient,  water  generated  in  the  catalyst 
layer  tends  to  remain  in  the  GDL  until  the  saturation  approaches 
0.2.  Owing  to  the  low  saturation  condition  in  the  gas  channel, 
the  saturation  in  the  GDL  hardly  exceeds  0.2  because  the  liq¬ 
uid  leaves  the  GDL  readily  once  this  threshold  value  is  attained. 


Fig.  4.  Capillary  diffusion  coefficient  used  in  the  calculation. 


Fig.  5.  Comparison  of  experimental  data  and  numerical  predictions  using  the 
user  Dcap  correlation  for  capillary  diffusion  coefficient. 

Therefore  the  second  half  of  the  correlation  for  capillary  diffu¬ 
sion  coefficient  (for  s  >  0.2)  does  not  play  a  role  in  the  numerical 
prediction. 

Fig.  7(a)  and  (b)  shows  the  current  density  profiles  of  experi¬ 
mental  data  and  numerical  predictions  using  different  correlation 
of  capillary  diffusion  coefficient.  The  discrepancy  in  the  cur¬ 
rent  density  prediction  is  obvious  for  low  humidification  cases 
(RH  25%  and  RH  0%  cases).  One  can  see  for  the  cases  tested, 
that  the  capillary  diffusion  coefficient,  or  rather  the  transport  of 
liquid  water,  has  little  impact  on  current  density  predictions. 
It  should  be  noted  that  the  sensitivity  to  RH  in  the  experi¬ 
ments  is  mostly  confined  to  the  inlet/outlet  regions  and  this 
is  not  entirely  surprising  as  such  counterflow  anode/cathode 
arrangements  can  promote  humidification  via  internal  water 
recirculation  as  described  by  Biichi  and  Srinivasan  [31]. 

3.4.  Membrane  properties 

3.4.1.  Parametric  study:  EOD 

In  order  to  assess  the  impact  of  the  electro-osmotic  drag 
coefficient,  a  parametric  study  is  undertaken,  and  the  results 
are  shown  in  Fig.  8.  The  diffusion  coefficient  for  water  in  the 
membrane  was  maintained  constant  (1  .2x10  12m2s  l),  while 
the  drag  coefficient,  n&  =  (n^/ 22) X,  was  varied  for  n £  from  a 
value  of  2.5  given  by  Springer  et  al.  [32]  to  one  order  mag¬ 
nitude  smaller.  It  should  be  noted  that  a  n%  value  of  unity 
was  reported  by  Zawodzinski  et  al.  [33]  and  used  in  simula¬ 
tion  by  some  researchers,  e.g.  Berg  et  al.  [13].  As  the  EOD  is 
decreased,  less  water  is  dragged  from  anode  to  cathode  near 
the  cathode  inlet,  whereas  near  the  anode  inlet  more  water  dif¬ 
fuses  across  the  membrane  from  cathode  to  anode.  This  results 
in  better  membrane  humidification  near  the  cathode  inlet,  and 
hence  a  higher  current  density,  see  Fig.  9.  A  lower  current  den¬ 
sity  region  results  near  the  anode  inlet,  and  overall  reducing  the 
EOD  yields  predicted  current  density  profiles  that  deviate  more 
from  experimental  data. 

3.4.2.  Sorption  isotherm 

The  predicted  water  mass  in  the  ME  A  obtained  with  var¬ 
ious  models  is  compared  to  measurements  in  Fig.  10.  The 
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Saturation 


(b)  Case  1  (RH100),  user  Dcap  capillary  diffusion  coefficient 


Saturation 

0.17 


(d)  Case  1 1  (RHO),  user  Dcap  capillary  diffusion  coefficient 

Fig.  6.  Predicted  saturation  in  the  cathode  side  of  the  MEA  and  gas  channel,  (a)  Case  1  (RH100),  default  capillary  diffusion  coefficient,  (b)  Case  1  (RH100),  user 
Dcap  capillary  diffusion  coefficient,  (c)  Case  1 1  (RHO),  default  capillary  diffusion  coefficient,  (d)  Case  1 1  (RHO),  user  Dcap  capillary  diffusion  coefficient. 


predicted  water  mass  plotted  with  triangles  is  calculated  using 
the  original  sorption  isotherm  given  in  Springer  et  al.  [32];  the 
square  symbols  represent  predictions  using  a  modified  sorption 
isotherm  that  has  a  smooth  transition  over  X  =  14-22  to  near 
unity  water  activity  to  account  for  the  different  water  content 
value  between  vapor  and  liquid  water  equilibrated  states  (i.e.  the 
so-called  Schroeder’s  paradox).  The  predicted  mass  of  water  in 
the  membrane  phase  is  essentially  constant  except  for  a  drop 
near  the  anode  inlet  due  to  the  low  RH  in  the  anode  channel. 
The  amount  of  liquid  water,  resulting  mainly  from  water  con¬ 
densation  in  the  cathode  side,  increases  along  the  channel  and 
slightly  near  the  anode  inlet  due  to  back  diffusion  as  discussed 
earlier.  The  experimental  data  falls  between  the  predictions  using 
the  Springer  isotherm  and  the  modified  isotherm.  Fig.  1 1  shows 
a  comparison  of  predicted  current  density  profiles  using  dif¬ 
ferent  sorption  isotherms  with  experimental  data.  Unlike  the 
case  for  MEA  water  mass,  the  current  density  prediction  using 
Springer’s  isotherm  yields  a  closer  match  to  the  experimental 
data. 


From  the  exercise  with  the  capillary  diffusion  coefficient,  it 
is  found  that  the  computed  water  content  in  the  inlet  region  is 
always  lower  than  measured  data.  It  is  then  reasonable  to  con¬ 
sider  possible  non-equilibrium  effects  in  the  membrane  at  low 
humidity  conditions.  Sorption  isotherms  are  obtained  by  equi¬ 
librating  the  membrane  over  a  sufficiently  long  period  of  time 
with  water  vapor  of  known  RH.  It  is  speculated  that  an  equilib¬ 
rium  state  might  not  be  reached  for  low  humidity  conditions  in 
an  operating  fuel  cell.  Therefore,  a  heuristic  approach  is  taken 
to  illustrate  how  non-equilibrium  effects  can  impact  the  predic¬ 
tions  of  water  content  and  current  density  profiles.  Fig.  12  shows 
the  default  sorption  isotherm  for  Nation  used  in  the  calculation 
versus  an  hypothetical  correlation  having  a  parabolic  form. 

Fig.  13  shows  predicted  water  content  profiles  using  different 
sorption  isotherms.  For  the  “non-equilibrium”  isotherm,  one  can 
see  an  increase  in  water  content  prediction,  particularly  in  the 
inlet  region.  This,  however,  shows  a  gap  between  the  numeri¬ 
cal  prediction  and  data.  If  one  shifts  the  prediction  curve  up  to 
make  the  plateau  region  of  the  curve  to  be  5  mg  cm-2,  the  inlet 
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Fig.  7.  Comparison  of  predicted  current  density  profiles  with  experimental  data 
(a)  CFD  versus  experimental  data  (b)  Predicted  current  density  profile  using 
default  and  userDcap  (marked  as  UserDcap  in  legend)  correlation  for  capillary 
diffusion  coefficient. 

region  of  the  water  content  compares  closely  to  experimental 
data,  whereas  the  cathode  outlet  portion  of  the  curve  shows  a 
large  discrepancy.  The  drop-off  of  water  content  in  cathode  out¬ 
let  shown  in  experimental  data  may  be  simulated  by  lowering 
the  coolant  flow  rate  as  will  be  discussed  later.  The  use  of  the 
“non-equilibrium”  sorption  isotherm  yields  a  more  satisfactory 


Fig.  9.  Comparison  current  density  data  with  numerical  predictions  for  different 
EOD  coefficient  values. 


Fig.  10.  Comparison  of  data  of  water  mass  in  the  MEA  with  numerical  predic¬ 
tions. 

result  in  current  density  prediction  than  in  water  content  profile, 
cf.  Fig.  13. 

3.4.3.  Water  dijfusivity  in  membrane 

For  the  sake  of  simplicity,  the  water  diffusivity  in  the  mem¬ 
brane  is  fixed  at  1.2E— 10  m2  s-1  for  most  of  the  calculations  in 
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Fig.  12.  Sorption  isotherms  used  in  calculation;  Nafion  =  default,  artifi¬ 
cial  =  hypothetical  correlation. 

the  present  study,  which  differs  from  that  given  in  Springer  et  al. 
[32].  This  value  is  indeed  close  to  the  diffusivity  corresponding 
to  an  activity  of  unity  in  the  original  Springer  paper.  Fig.  14(a) 
and  (b)  shows  that  the  difference  in  predicted  water  content  and 
current  density  using  the  constant  diffusivity  and  the  correlation 
of  Springer  et  al.  is  not  significant. 

3.5.  Catalyst  layer 

3.5.1.  Specific  surface  area  for  fitting  stoichiometry 
sensitivity  data 

The  stoichiometry  sensitivity  of  a  unit  cell  is  important  for 
fuel  cell  stack  operation  because  in  a  stack  the  cells  are  likely 
to  operate  at  different  stoichiometric  ratios  and  minimum  loss 
in  performance  due  to  unevenly  distributed  stoichiometric  ratio 
is  desired.  Therefore,  in  addition  to  the  single  cell  performance, 
the  stoichiometry  sensitivity  is  an  important  measure  of  a  unit 
cell.  As  the  stoichiometric  ratio  approaches  unity,  performance 
of  the  cell  decreases  due  to  mass  transfer  limitation,  which  can 
occur  at  different  region  of  the  unit  cell,  e.g.  in  the  GDL  due  to 
liquid  water  blockage  or  in  the  catalyst  layer  due  to  pore  level 
mass  diffusion.  The  effects  due  to  mass  transfer  through  the  GDL 


Fig.  13.  Water  content  predictions  using  different  sorption  isotherms  and  aver¬ 
age  pore  size. 


Water  content  of  ME  A 


Puck  # 

Current  density  profile 

Fig.  14.  Comparison  of  water  content  and  current  density  predictions  using  a 
constant  water  diffusivity  in  membrane  and  the  Springer  model:  (a)  water  content 
of  ME  A  and  (b)  current  density  profile. 

manifest  in  the  saturation  transport,  which  has  been  included  in 
the  aforementioned  tests.  The  pore-level  mass  transfer  in  the 
catalyst  layer  is  treated  via  the  so-called  “diffusion-reaction  bal¬ 
ance,”  cf.  [21,27].  The  model  parameter  used  to  fit  experimental 
data  of  stoichiometric  sensitivity  is  the  specific  surface  area  (S/V) 
of  the  catalyst.  Fig.  15  shows  a  comparison  of  data  of  cell  volt¬ 
age  with  numerical  predictions  using  S/V  ranging  from  2  to  100. 
It  is  found  that  the  value  S/V  =  3  fits  the  data  best. 

3.5.2.  On  the  limiting  current  predicted  using  the  current 
CFD  model 

The  limiting  current  for  a  unit  cell  is  dependent  on  many  fac¬ 
tors  along  the  pathway  of  oxidant  transfer.  From  the  gas  channel 
inlet  to  the  catalyst  surface,  transfer  of  oxygen  is  subjected  to 
limiting  mechanisms  including  (a)  oxygen  concentration  gradi¬ 
ent  along  the  channel  due  to  local  consumption,  (b)  convective 
mass  transfer  from  the  channel  bulk  flow  to  the  GDL  surface, 
(c)  diffusive  transfer  through  the  GDL  in  the  presence  of  porous 
solid  and  liquid  water,  (d)  dissolution  of  oxygen  into  water  and 
electrolyte,  (e)  diffusion  of  oxygen  in  the  and  mass  diffusion 
in  the  pore  level,  and  (f)  diffusion  of  dissolved  oxygen  in  the 
electrolyte.  In  the  current  CFD  model  for  PEMFC,  factors  (a), 
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Fig.  15.  Predicted  cell  voltage  as  a  function  of  oxygen  stoichiometric  ratio  with 
various  S/V  (sov  in  the  legend)  values. 

(b),  (c)  are  included  in  the  simulation.  Contribution  by  factors 
(d)  and  (f)  are  not  explicitly  modeled,  however,  both  factors  can 
be  lumped  into  the  reaction-diffusion  balance  (e),  and  can  be 
handled  by  one  control  parameter,  e.g.  the  specific  surface  area 
of  catalyst  S/V.  Ultimately  the  transport  processes  that  have  not 
been  accounted  for  should  be  modeled  as  more  experimental 
data  and  analysis  become  available. 

To  gain  insight  into  the  limiting  current  change  due  to  the 
model  parameter  S/V,  a  series  of  polarization  curve  predictions 
are  calculated  for  various  stoichiometric  factors  for  two  oxygen 
concentrations.  Fig.  16  shows  a  comparison  of  the  numerical 
results.  The  stoichiometric  factor  is  calculated  based  on  air  at 
I-  1  A  cm-2.  The  stoichiometric  ratio  {stoich  in  Fig.  16)  of  these 
cases  is  adjusted  by  changing  the  inlet  gas  velocity;  e.g.  for  sto¬ 
ich  =  1 . 1 ,  the  inlet  gas  velocity  is  1 . 1  times  the  inlet  velocity 
required  for  an  average  current  density  of  1  A  cm-2  when  air  is 
used.  The  cases  with  high  oxygen  concentration  has  the  oxygen 
and  nitrogen  concentration  swapped  for  air.  For  the  air  cases, 


Stoich  defined  as  flow  rate  versus  1=1  A/cm2  air 


Fig.  16.  Predicted  polarization  curves  with  various  oxygen  flow  rates  and 
concentrations.  The  oxygen  stoichiometric  factor  is  based  on  using  air  at 
1=  1  A  cm-2. 


DP  =  Dew  Point 


Fig.  17.  Numerical  predictions  of  mass  flow  rate  for  different  anode  humidifi¬ 
cation. 


the  limiting  current  when  a  high  S/V  is  used  coincides  with  the 
stoichiometry  of  oxygen.  When  S/V  =  3  is  used,  which  provides 
the  best  fit  for  stoichiometry  sensitivity  data,  the  limiting  cur¬ 
rent  shifts  to  a  lower  value  than  the  high  S/V  case.  For  the  higher 
oxygen  concentration  cases,  zero  cell  potential  occurs  prior  to 
the  limiting  current  condition.  It  should  be  noted  that  for  the 
cases  with  high  oxygen  concentration  ( Xo2  =  0.73  in  Fig.  16), 
the  amount  of  oxygen  flow  rate  at  inlet  is  in  fact  equivalent  to 
stoich  =  3.46  based  on  air  because  of  the  high  oxygen  concen¬ 
tration  at  the  inlet.  The  implication  of  this  study  is  that  one  may 
use  limiting  current  data  to  calibrate  the  model  parameter  used 
for  the  reaction-diffusion  balance. 

3.6.  Operating  conditions 

3.6.1.  Anode  humidification 

Fig.  17  shows  the  flow  rates  for  four  different  anode  humidifi¬ 
cation  conditions  corresponding  to  inlet  dew  point  temperatures 
of  58,  63,  70,  and  75  °C.  The  only  noticeable  difference  is  in  the 
anode  water  vapor  flow,  which  exhibits  a  monotonic  decrease 
along  the  channel  except  for  the  drier  case  (58  °C  dew  point).  It 
is  interesting  to  note  that  near  the  anode  outlet  the  water  vapor 
flow  rates  all  converge  to  a  similar  level.  The  cases  shown  in 
Fig.  18  are  calculated  with  the  same  potential  difference  across 
the  unit  cell.  The  effect  due  to  anode  humidification  is  there- 


Fig.  18.  Numerical  predictions  of  current  density  and  relative  humidity  for 
different  anode  humidification. 
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Fig.  19.  Predicted  coolant  temperature  profiles  in  the  axial  direction  for  various 
coolant  flow  rates. 


Coolant  flow  rate  (slpm) 

Fig.  21.  Temperature  difference  between  coolant  inlet  and  outlet  at  various 
coolant  flow  rates. 


fore  reflected  in  the  change  of  local  current  density  near  the 
anode  inlet  as  shown  in  Fig.  18,  which  increases  with  anode 
humidification. 

3.6.2.  Dependence  of  current  density  and  water  content 
predictions  on  coolant  flow  rate 

In  comparing  the  simulation  results  with  experimental  data, 
the  coolant  flow  rate  needs  to  be  adjusted  in  order  to  obtain  tem¬ 
perature  gradient  in  the  axial  direction  similar  to  that  recorded 
in  the  experiments  because  of  uncertainty  issue  in  actual  exper¬ 
iment  operation.  Fig.  19  shows  the  temperature  profiles  in  the 
axial  direction  as  a  function  of  inlet  coolant  flow  rate.  As  the 
coolant  flow  rate  is  reduced,  the  temperature  profiles  rise  and 
temperature  gradient  between  coolant  inlet  and  outlet  becomes 
higher.  The  temperature  difference  between  coolant  inlet  and 
outlet  is  shown  in  Fig.  20.  The  sensitivity  of  the  cell  perfor¬ 
mance  appears  to  change  drastically  when  the  coolant  flow  rate 
is  reduced  to  certain  point,  cf.  ca.  0. 1  slpm.  Similar  to  stoichiom¬ 
etry  sensitivity,  the  cell  performance  is  also  sensitive  to  coolant 
flow  rate,  or  rather  the  heat  removal  capacity  defined  as  mCp. 
As  is  show  in  Fig.  21,  once  the  heat  removal  capacity  drops 
below  certain  threshold  value,  the  cell  performance  will  be  sig¬ 
nificantly  reduced.  Such  condition  may  occur  when  good  flow 


Fig.  20.  Predicted  water  content  and  current  density  profiles  for  various  coolant 
flow  rates. 


0.58 

0.56 

r 

0.54 

1 

I 

> 

0.52 

1 

1 

1 

1 

> 

0.5 

1 

1 

0.48 

1 

_  1 

1 

-  « 

0.46 

0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8 

Coolant  flow  rate  (slpm) 

Fig.  22.  Predicted  cell  voltage  as  a  function  of  coolant  flow  rate. 

sharing  of  the  coolant  in  a  stack  is  not  maintained.  Fig.  22  shows 
the  predicted  water  profiles  and  current  density  profiles  at  vari¬ 
ous  coolant  flow  rates.  When  the  coolant  flow  rate  is  low,  the  high 
temperature  near  the  coolant  outlet  (located  on  the  same  side  as 
the  cathode  gas  outlet  in  the  MRED  configuration)  causes  mem¬ 
brane  dry-out  and  an  expansion  of  the  low  water  content  region 
near  the  outlet. 

4.  Conclusions 

The  CFD-based  computational  framework  described  in  Part 
1,  was  in  this  paper  systematically  validated  against  spatially 
resolved  measurements,  including  water  balance  and  current 
density  distributions,  obtained  under  a  broad  range  of  operating 
conditions  that  stretch  the  computational  model  in  general,  and 
the  embedded  constitutive  relations  of  some  of  the  sub-models 
in  particular.  While  simulation  results  based  on  membrane  prop¬ 
erties  reported  by  Springer  et  al.  [32]  capture  the  overall  trends, 
examination  of  the  predicted  water  mass  in  the  MEA  using 
various  formulations,  show  that  experimental  data  is  bracketed 
by  predictions  using  the  “standard”  sorption  isotherm  correla¬ 
tion  and  predictions  using  a  modified  correlation  that  takes  into 
consideration  liquid  equilibrated  conditions. 
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Based  on  the  analysis,  adjustments  were  proposed  for  some 
of  the  constitutive  equations,  including  the  capillary  diffusion 
coefficient  for  liquid  saturation  and  sorption  isotherms  for  water 
uptake  as  well  as  model  parameters  such  the  specific  surface 
area.  These  adjustments  were  quite  effective  in  improving  the 
overall  fit  between  experiments  and  simulations.  In  the  limit  of 
low  RH  conditions,  it  is  found  that  adjustments  in  the  coolant 
flow  rate  are  required  to  achieve  a  good  fit  with  experimental 
with  respect  to  water  content  in  the  cathode  outlet  region.  In 
terms  of  water  transport  across  the  membrane,  the  cathode  out¬ 
let  region  (diffusion-dominated)  might  differ  from  the  cathode 
inlet  region  (drag-dominated)  in  the  direction  of  net  water  flux. 
It  is  nonetheless  possible  to  identify  some  critical  parameters 
that  yield  a  good  fit  of  the  water  content  profile  at  both  ends 
of  the  cell  simultaneously,  indicating  capabilities  to  account  for 
a  large  range  of  changes  in  terms  of  dominant  transport  mech¬ 
anisms  and  a  level  of  generality  that  is  promising  in  terms  of 
functionality  for  design.  The  extensive  validation  undertaken 
here  would  further  benefit  from  experimental  method  that  could 
accurately  differentiate  water  content  in  the  GDL  from  that  in 
the  membrane.  Finally,  it  is  clear  that  fundamental  studies  are 
required  to  set  some  of  the  constitutive  relations  and  correlations 
on  firmer  ground.  Some  very  recent  developments  for  charac¬ 
terizing  the  two-phase  flow  in  gas  diffusion  media  [29,30,34] 
have  resulted  for  instance  in  new  capillary  and  relative  perme¬ 
ability  functions  that  should  further  enhance  the  reliability  and 
generality  of  computational  fuel  cell  models. 
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Appendix  A.  Sample  calculation  for  the  mass  of  water 
in  the  MEA 

For  a  GDL  of  250  jxm  in  thickness  and  porosity  0.6,  and  a 
membrane  of  50  |xm. 

•  For  a  saturation  of  0.2  in  the  GDL,  the  water  content 
is  250E— 6  m  x  0.6  x  0.2  x  1000  kg  m-3  x  lE6mg  kg-1  x 
IE— 4  m2  cm  2  =  3  mg  cm-2. 

•  For  a  fully  humidified  membrane,  the  water  content  is 
A,  x  50E— 6mx  1980kgm_3-l.l  kg  mol-1  x  0.018kg 
mol-1  x  lE6mg  kg-1  x  IE— 4  m2  cm'2  =0.1  62a  mg 
cm  z. 


•  For  A  =  14  (vapor-equilibrated),  the  value  of  water  content 
is  about  2.3  mg  cm-2  and  3.6  mg  cm-2  for  A.  =  22  (liquid- 
equilibrated). 
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